using System;
using System.Collections;
using System.Collections.Generic;
using System.Text;


namespace LateralProcesses
{
    /// <summary>
    /// This is a HBC-EC CDU with lateral flow routing between CDUs
    /// </summary>
    public class CDUHBVEC: CDUClass
    {

        #region private variables

        // Compartment members
        private Reservoir _fastReservoir = new Reservoir();
        private Reservoir _slowReservoir = new Reservoir();   
        
        // private variables
        private double _evaporation;            // evaporation
        private double _testV1;                 // test a variable for debugging
        private double _testV2;                 // test a variable for debugging

        // model parameter members
        private string routingConfig = "paralell";
        private double KF; //fast reservoir coefficient (1/d)
        private double ALPHA; //fast reservoir exponent
        private double KS; //slow reservoir coefficient (1/d)
        private double FRACPERC; // fraction of runoff directed to fast reservoir or percolation in serial configuration

        private double _initFastReservoirDischarge;
        private double _initSlowReservoirDischarge;
        private double _landArea; //m2
        private double _lakeArea; //m2
        private double _glacierArea; //m2
        private double _glacierRunoff;
        private double _lakeRunoff;
        private double timestepSeconds;
        private double timeStepHrs;

        #endregion


        #region CDUGeorg Public Methods

        public override void Initialize(Hashtable InitialConditions)
        {
            throw new Exception("The method or operation is not implemented.");
            //initialize(InitialConditions);
        }

        //here is the actual modeling
        //thats where all the hard work is done.
        public override void Update(DateTime TimeStep, Dictionary<string, double> GruToCdu)
        {
            double landRunoff = CalcWaterInputRate(GruToCdu);
            double fastResElev = _fastReservoir.RElevation;
            double fastResDelta = 0.0;
            double slowResElev = _slowReservoir.RElevation;
            double slowResDelta = 0.0;

            if (routingConfig.ToUpper() == "PARALELL")
            {
                //fast response reservoir
                double startInflow = _fastReservoir.StartInflow; //[mm]
                double endInflow = FRACPERC * landRunoff + lateralInflowFast / Area;//add inflow [mm] from parent cdu
                double initialFastElev = _fastReservoir.RElevation;
                fastResElev = RoutingFastRK3(startInflow, endInflow, initialFastElev);
                
                // update state variables
                _fastReservoir.StartInflow = endInflow;

                //slow response reservoir
                startInflow = _slowReservoir.StartInflow;
                endInflow = (1 - FRACPERC) * landRunoff + lateralInflowSlow / Area; //add inflow [mm] from parent cdu
                double initialSlowElev = _slowReservoir.RElevation;
                slowResElev = RoutingSlowRK3(startInflow, endInflow, initialSlowElev);
                
                // update state variables
                _slowReservoir.StartInflow = endInflow;
                
            }

            if (routingConfig.ToUpper() == "SERIAL")
            {

                //fast response reservoir
                double startInflow = _fastReservoir.StartInflow;
                double endInflow = landRunoff + lateralInflowFast / Area; //add inflow [mm] from parent cdu
                double initialFastElev = _fastReservoir.RElevation;
                fastResElev = RoutingFastRK3(startInflow, endInflow, initialFastElev);

                // update state variables
                _fastReservoir.StartInflow = endInflow;

                //Calculate percolation from fast to slow reservoir
                double perc = GetProratedParameter(FRACPERC);
                if (fastResElev > perc)
                    fastResElev -= perc;
                else
                {
                    perc = fastResElev;
                    fastResElev = 0.0;
                }

                //slow response reservoir
                startInflow = _slowReservoir.StartInflow;
                endInflow = perc + lateralInflowSlow / Area;  //add inflow [mm] from parent cdu;
                double initialSlowElev = _slowReservoir.RElevation;
                slowResElev = RoutingSlowRK3(startInflow, endInflow, initialSlowElev);
                
                //set current percolation (start percolation for next time step)
                _slowReservoir.StartInflow = endInflow;
                     
            }

             //storage elevations
            _fastReservoir.RElevation = fastResElev;
            _slowReservoir.RElevation = slowResElev;
            
            //runoff [l/timestep]
            _fastReservoir.RRunoff = ComputeFastReservoirRunoff(fastResElev) * Area;
            _slowReservoir.RRunoff = ComputeSlowReservoirRunoff(slowResElev) * Area;

          

            // output [l/timestep]
            RunoffFast = _fastReservoir.RRunoff;
            RunoffSlow = _slowReservoir.RRunoff;
            RunoffTotal = _fastReservoir.RRunoff + _slowReservoir.RRunoff;
            _testV1 = fastResElev;
           
            lateralInflowTotal = lateralInflowFast / Area + lateralInflowSlow / Area;
            _testV2 = slowResElev;
            //reset lateral inflows
            lateralInflowFast = 0.0;
            lateralInflowSlow = 0.0;
            updateStateAndFlux();
        }

  
        /// <summary>
        /// Route the fast reservoir using a 3rd order Runga-Kutta method
        /// </summary>
        /// <param name="inflowStart"></param>
        /// <param name="inflowEnd"></param>
        /// <param name="resElev"></param>
        /// <returns></returns>
        private double RoutingFastRK3(double inflowStart, double inflowEnd, double resElev)
        {

            double I1 = inflowStart;
            double I2 = inflowStart + (inflowEnd - inflowStart) / 3.0;
            double I3 = inflowStart + 2.0 * (inflowEnd - inflowStart) / 3.0;

            double H1 = I1 - ComputeFastReservoirRunoff(resElev);
            double H2 = I2 - ComputeFastReservoirRunoff(resElev + H1 / 3.0);
            double H3 = I3 - ComputeFastReservoirRunoff(resElev + 2.0 * H2 / 3.0);

            double deltaH = (0.25 * H1) + (0.75 * H3);

            return resElev + deltaH;
        }
        /// <summary>
        /// Route the slow reservoir using a 3rd order Runga-Kutta method
        /// </summary>
        /// <param name="inflowStart"></param>
        /// <param name="inflowEnd"></param>
        /// <param name="resElev"></param>
        /// <returns></returns>
        private double RoutingSlowRK3(double inflowStart, double inflowEnd, double resElev)
        {

            double I1 = inflowStart;
            double I2 = inflowStart + (inflowEnd - inflowStart) / 3.0;
            double I3 = inflowStart + 2.0 * (inflowEnd - inflowStart) / 3.0;

            double H1 = I1 - ComputeSlowReservoirRunoff(resElev);
            double H2 = I2 - ComputeSlowReservoirRunoff(resElev + H1 / 3.0);
            double H3 = I3 - ComputeSlowReservoirRunoff(resElev + 2.0 * H2 / 3.0);

            double deltaH = (0.25 * H1) + (0.75 * H3);

            return resElev + deltaH;
        }
        #endregion


        /// <summary>
        /// This is an example of a derived class from (CDUClass) which computes a
        /// one-dimensional water balance for a CDU. Note that we must override all 
        /// abstract methods.
        /// </summary>
        public CDUHBVEC(double timeStepHrs, string id, string idDown, List<string> cduId, List<double> cduArea, double area,
            double Z, double slope, double aspect, string routingConfig, bool output, Dictionary<string,double> cduGlobalParameters, Dictionary<string,double> cduLocalParameters)
            : base(id, idDown, cduId, cduArea, area, Z, slope, aspect, output, cduGlobalParameters, cduLocalParameters)
        {
            this.timeStepHrs = timeStepHrs;
            this.timestepSeconds = timeStepHrs * 60 * 60;
            this.routingConfig = routingConfig;
            assignParameters();

            // initialize state variables (later add initial condidtions)
            RunoffFast = 0.0;
            _evaporation = 0.0;
            RunoffFast = 0.0;
            RunoffSlow = 0.0;
            lateralInflowFast = 0.0;
            lateralInflowSlow = 0.0;

            // initialize Reservoirs (later on move to initialize())
            _fastReservoir.RName = "FastReservoir";
            _slowReservoir.RName = "SlowReservoir";

            // Set initial reservoir elevations from discharge
            _fastReservoir.RElevation = ComputeFastReservoirElevation(_initFastReservoirDischarge);
            _slowReservoir.RElevation = ComputeSlowReservoirElevation(_initSlowReservoirDischarge);
            updateStateAndFlux();
        }

          
        #region Private methods


        /// <summary>
        /// This function retrieves the parameters from the Parameters hashtable and assigns them to 
        /// 'real' variables
        /// </summary>
        private void assignParameters()
        {
            KF = base.CduGlobalParameters["kf"];
            ALPHA = base.CduGlobalParameters["alpha"];
            KS = base.CduGlobalParameters["ks"];
            FRACPERC = base.CduGlobalParameters["fracperc"];
            _initFastReservoirDischarge = base.CduLocalParameters["inifastq"];
            _initSlowReservoirDischarge = base.CduLocalParameters["inislowq"];
        }

        private void updateStateAndFlux()
        {
            FluxStateVars = new Dictionary<string, double>();
            FluxStateVars.Add("runofftotal", RunoffTotal / timestepSeconds);
            FluxStateVars.Add("inflowtotal", lateralInflowTotal);
            FluxStateVars.Add("runofffast", RunoffFast / timestepSeconds);
            FluxStateVars.Add("inflowfast", lateralInflowFast);
            FluxStateVars.Add("runoffslow", RunoffSlow / timestepSeconds);
            FluxStateVars.Add("inflowslow", lateralInflowSlow);
            FluxStateVars.Add("evaporation", _evaporation);
            FluxStateVars.Add("testV1", _testV1);
            FluxStateVars.Add("testV2", _testV2);
        }
        /// <summary>
        /// calculates a reservoir elevation to achieve a given 
        /// reservoir outflow m3/timestep
        /// </summary>
        /// <param name="outflow"></param>
        /// <returns>
        /// fast reservoir elevation in millimeters
        /// </returns>
        private double ComputeFastReservoirElevation(double outflow)
        {
            double KFPR = GetProratedParameter(KF);
            return Math.Pow(outflow / Area * 1000.0 / KFPR, 1 / (1 + ALPHA)); 
        }

        private double ComputeSlowReservoirElevation(double outflow)
        {
            double KSPR = GetProratedParameter(KS);
            return outflow / Area * 1000.0 / KSPR;
        }
     
        /// <summary>
        /// Computes the water input rate for CDU as 
        /// area weigthed mean of GRUs within the CDU
        /// </summary>
        private double CalcWaterInputRate(Dictionary<string, double> GruToCdu)
        {
            double waterInputRate = 0.0;
            for (int i = 0; i < _gruId.Count; i++)
            {
                waterInputRate += GruToCdu[_gruId[i]] * _gruArea[i];
                
            }

            return waterInputRate / Area;
        }
        /// <summary>
        /// Calculates the amount of runoff from fast reservoir 
        /// </summary>
        /// <param name="elevation">elevation of the reservoir [mms]</param>
        /// <returns>reservoir runoff in millimetres</returns>
        private double ComputeFastReservoirRunoff(double elevation)
        {
            
            double runoff = GetProratedParameter(KF) * Math.Pow(elevation, 1 + ALPHA);

            // high values of alpha may actually withdraw more water from reservoir than
            // available. Take only what is avaialble 
            return Math.Min(runoff,elevation);

        }
        /// <summary>
        /// Calculates the amount of runoff from slow reservoir, 
        /// in mms given the elevation of the reservoir, in mms 
        /// Returns reservoir runoff in millimetres
        /// </summary>
        /// <param name="elevation"></param>
        /// <returns></returns>

        private double ComputeSlowReservoirRunoff(double elevation)
        {
            return GetProratedParameter(KS) * elevation;
        }

        /// <summary>
        /// Returns Parameter (1/day) prorated to the timestep
        /// </summary>
        /// <param name="value"></param>
        /// <returns></returns>
        private double GetProratedParameter(double value)
        {
            if (timestepSeconds != 86400)
            {
                return value / 86400 * timestepSeconds;
            }
            return value;
        }
        #endregion


    }
}
